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Abstract —  An  accurate,  computationally  efficient  and 
fully-automated  algorithm  for  the  alignment  of  2D  seri¬ 
ally  acquired  sections  forming  a  3D  volume  is  presented. 
The  method  accounts  for  the  main  shortcomings  of  3D  im¬ 
age  alignment:  corrupted  data  (cuts  and  tears),  dissimilar¬ 
ities  or  discontinuities  between  slices,  non  parallel  or  miss¬ 
ing  slices.  The  approach  relies  on  the  optimization  of  a 
global  energy  function,  based  on  the  object  shape,  measur¬ 
ing  the  similarity  between  a  slice  and  its  neighborhood  in 
the  3D  volume.  Slice  similarity  is  computed  using  the  dis¬ 
tance  transform  measure  in  both  directions.  No  particular 
direction  is  privileged  in  the  method  avoiding  global  off¬ 
sets,  biases  in  the  estimation  and  error  propagation.  The 
method  was  evaluated  on  real  images  (medical  and  biologi¬ 
cal  3D  data)  and  the  experimental  results  demonstrated  the 
method’s  accuracy  as  reconstuction  errors  are  less  than  1 
degree  in  rotation  and  less  than  1  pixel  in  translation. 

Keywords:  serially  acquired  images,  misalignment,  image 
registration,  registration  error,  non-overlapping  structures, 
pixel  similarity  measure,  deterministic  optimization. 

I.  Introduction 

Three-dimensional  reconstruction  of  medical  images  (tis¬ 
sue  sections,  CT  and  autoradiographic  slices)  is  now  an  in¬ 
tegral  part  of  biomedical  research.  Reconstruction  of  such 
data  sets  into  3D  volumes,  via  the  registrations  of  2D  sec¬ 
tions,  has  gained  an  increasing  interest.  The  registration 
of  multiple  slices  is  of  utmost  importance  for  the  correct 
3D  visualization  and  morphometric  analysis  (e.g.  surface 
and  volume  representation)  of  the  structures  of  interest. 
Several  alignment  algorithms  have  been  proposed  in  that 
framework.  A  review  of  general  medical  image  registration 
methods  is  presented  in  [1],  [2],  [3]. 

The  principal  3D  alignment  (reconstruction  from  2D  im¬ 
ages)  methods  may  be  classified  in  the  following  categories: 
fiducial  marker-based  methods  [4],  feature-based  methods 
using  contours,  crest  lines  or  characteristic  points  extracted 
from  the  images  [5],  [6],  and  gray  level-based  registration 
techniques  using  the  intensities  of  the  whole  image  [7] ,  [8] , 
[9],  [10].  Most  of  the  above  mentioned  techniques  do  not 
simultaneously  consider  the  two  major  difficulties  involved 
in  medical  and  CT  scanned  data  registration. 

At  first,  consecutive  slices  may  differ  significantly  due  to 
distortions,  discontinuities  in  anatomical  structures,  cuts 
and  tears.  These  effects  are  more  pronounced  when  distant 
slices  are  involved  in  the  registration.  From  this  point  of 
view,  a  registration  method  must  be  robust  to  missing  data 


or  outliers  [7],  [10]. 

Besides,  registering  the  slices  sequentially  (the  second 
with  respect  to  the  first,  the  third  with  respect  to  the  sec¬ 
ond,  etc.)  leads  to  different  types  of  misregistration.  If  an 
error  occurs  in  the  registration  of  a  slice  with  respect  to 
the  preceding  slice,  this  error  will  propagate  through  the 
whole  volume.  Also,  if  the  number  of  slices  to  be  registered 
is  large,  a  global  offset  of  the  volume  may  be  observed,  due 
to  error  accumulation  [8]. 

In  this  paper,  a  solution  to  the  above  mentioned  short¬ 
comings  is  presented.  A  global  energy  function  having  as 
variables  the  rigid  transformation  parameters  (2D  transla¬ 
tion  and  rotation)  of  a  given  slice  with  respect  to  a  local 
symmetric  neighborhood  is  proposed.  Global  energy  func¬ 
tions  are  a  powerful  tool  in  computer  vision  applications 
but  they  have  not  yet  been  considered  for  the  registration 
of  serially  acquired  slices. 

Our  approach  was  inspired  by  the  technique  presented  in 
[11],  which  consists  in  minimizing  a  global  energy  function 
with  the  Iterative  Closest  Point  algorithm  [12],  to  regis¬ 
ter  multiple,  partially  overlapping  views  of  a  3D  structure. 
The  global  energy  function  implemented  in  our  approach 
is  associated  with  a  pixel  similarity  metric  based  on  the 
Euclidean  distance  transform  [13]. 

The  remainder  of  the  paper  is  organized  as  follows.  The 
global  energy  function  formulation  and  the  associated  reg¬ 
istration  algorithm  is  presented  in  section  II,  experimen¬ 
tal  results  are  presented  in  section  III  and  conclusions  are 
drawn  in  section  IV. 

II.  A  GLOBAL  ENERGY  FUNCTION  FORMULATION 

Before  presenting  the  alignment  method,  the  notations 
used  in  our  formulation  are  introduced.  A  set  of  2D  serially 
acquired  slices  is  represented  by: 

V  =  {Ii\  i  =  1  ...N}  (1) 

where  Ii  is  a  slice  (a  2D  image)  and  N  denotes  the  total 
number  of  slices.  A  pixel  of  a  2D  slice  is  represented  by: 
p  =  (x,2/)t,  so  that  Ii{p)  corresponds  to  the  gray  level 
(intensity)  of  pixel  p  of  slice  i.  Nx  and  Ny  designate  the 
number  of  pixels  of  each  slice  in  the  horizontal  and  vertical 
direction  respectively. 

Standard  two-dimensional  rigid  alignment  consists  of  es¬ 
timating  the  rigid  transformation  parameters  (translation 
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tXj  ty  and  rotation  by  angle  0)  that  have  to  be  applied  to 
the  image  to  be  aligned  (floating  image)  in  order  to  match 
a  r ef#M<?<^M&ge? 3 rd  Annual  Conference  -  IEEE/EMBS  Oct.25 

In  the  approach  proposed  here,  the  alignment  of  the  2D 
sections,  within  the  3D  volume,  is  considered  globally  by 
minimizing  an  energy  function  E(-),  which  expresses  the 
similarity  between  the  2D  sections: 

N  N  Nx  x  Ny 

£(©)  =  EE  E  KhiTs^JATsM)))  (2) 

2=1  j  =  1  p=  1 

where  /(•)  is  a  similarity  metric,  R  denotes  slice  k  and  T®k 
designates  a  rigid  transformation  with  parameters  0^  = 
{tkX,ty,Ok}. 

Equation  (2)  indicates  that  for  a  given  set  of  rigid  trans¬ 
formation  parameters  T© . ,  applied  to  the  slice  to  be  aligned 
/i,  the  similarity  between  the  transformed  slice  R(T@i(p)) 
and  all  of  the  other  already  transformed  slices  Ij(T®j(p)) 
in  the  volume  is  accumulated  in  the  energy  function. 

Assuming  that  function  /(•)  is  symmetric: 

f (liiTe, (p)),  I j  (T&J  (p)))  =  f(Ij(Tej  (p)),  I i (To ^p)))  (3) 

which  is  the  case  for  the  pixel  similarity  functions  consid¬ 
ered  here,  yields  the  following  global  minimization  prob¬ 
lem: 

N  N  Nx  x  Ny 

Q  =  argmm'£Yl  E  fOi(TeM),  IjiTe^p)))  (4) 

2=1  j  =  1  p=  1 
j<i 

Without  additional  constrains,  the  optimization  problem 
(4)  has  clearly  an  infinite  number  of  solutions  (if  the  set  of 
rigid  transformations  ,  Tq2  ,  . . .  T^n  }  is  a  solution,  the 
same  holds  true  for  oTa,  X^oTa,  . . .  T^n  oTa},  where 

Ta  is  an  arbitrary  2D  rigid  transformation) .  To  remove  this 
ambiguity,  the  transformation  applied  to  an  arbitrary 
chosen  slice  k  is  constrained  to  the  identity  transformation 
(we  have  chosen  k  =  1  in  our  implementation).  As  a  result, 
there  are  3(N  —  1)  parameters  to  estimate. 

It  is  common  sense  that  distant  slices  present  very  little 
similarity  due  to  anatomy  and  it  would  be  more  appropri¬ 
ate  to  measure  the  energy  function  only  for  slices  presenting 
at  least  some  similarities.  Therefore,  the  support  region  of 
function  /(•)  has  been  limited  to  a  neighborhood  of  radius 
R  centered  at  each  slice  and  set: 

/(/^(p)), /^(p)))  =  0,  V  \i-j\>R  (5) 

Thus,  the  following  alignment  algorithm  is  associated 
with  the  energy  function  (4): 

•  do  until  convergence. 

—  declare  all  slices  unvisited. 

-  do  until  all  slices  are  declared  visited. 

*  randomly  chose  an  unvisited  slice  R  E  V. 

*  update  the  rigid  transformation  parameters  T©. 

bringing  into  alignment  slice  R  with  the  other 


slices  in  the  neighborhood  of  i,  by  minimization 
of  the  following  local  energy  function  (eq.  6). 

-28,  200|,  ed . 

—  end  do 

•  end  do 

N  N  NxxNy 

^(0*)d=E  E  E  fOiiTs^JATsM)))  (6) 

2=1  j  =  1  p=  1 

\i-j\<R 

The  minimization  of  the  local  energy  function  (4)  is  con¬ 
ducted  by  a  deterministic  optimization  algorithm,  known 
as  Iterated  Conditional  Modes  (ICM)  [14].  ICM  is  a  dis¬ 
crete  Gauss  Seidel-like  optimization  technique,  accepting 
only  configurations  decreasing  the  objective  function.  Let 
us  notice  that  the  parameter  0^  corresponding  to  the  min¬ 
imum  value  of  the  local  energy  function  Ei(QR  (Equ.  6) 
also  corresponds  to  a  local  minimum  value  of  the  global 
energy  function  E(Q)  with  respect  to  0^,  keeping  all  other 
parameters  0j,  j  ^  i  fixed.  It  is  thus  easy  to  see  that  the 
described  algorithm  converges  towards  a  local  minimum  of 
the  initial  energy  function  (2).  This  local  minimum  corre¬ 
sponds  to  a  satisfactory  registration,  since  the  initial  align¬ 
ment  of  the  2D  sections  is  generally  close  to  the  desired 
solution  (if  this  is  not  the  case,  a  good  initialization  may 
be  obtained  by  a  standard  coarse  alignment  technique  such 
as  principal  axes  registration).  It  is  thus  not  necessary  to 
resort  here  to  greedy  global  optimization  procedures,  such 
as  simulated  annealing  or  genetic  algorithms. 

Further  improvement  of  the  solution  is  obtained  by  a 
gradient  decent  technique.  To  speed  the  algorithm  up  a 
multigrid  data  processing  is  also  implemented. 

The  pixel  similarity  metric  associated  with  the  above  de¬ 
scribed  global  energy  function  is  based  on  a  distance  trans¬ 
form  ([13],  [15])  (also  known  as  chamfer  matching  technique 
[16])  and  it  is  computed  from  the  3D  object  contours  ([17]). 
A  distance  transformation  is  an  operation  that  converts  a 
binary  picture,  consisting  of  feature  and  non-feature  ele¬ 
ments  (contours) ,  to  a  picture  where  each  pixel  has  a  value 
that  approximates  its  distance  to  the  nearest  contour  point. 

Thus,  using  the  distance  transform  D(p)  of  the  reference 
slice  the  method  aligns  the  floating  slice  by  minimizing  the 
distance  between  the  contours  of  the  images.  For  further 
details  of  the  chamfer  matching  method  the  reader  may 
refer  to  [16]. 

Considering  the  slices  per  triplets,  which  is  very  common 
for  standard  reconstruction  problems  (i.e.  setting  R=  1  in 
eq.  5),  the  estimation  of  the  alignment  parameters  0  in¬ 
volves  the  non-linear  similarity  metric: 

f(TeM)  =  A-i(T0i_1(p))  +  A+i(Tei+1(p)), 

UiTeM)^  0 

where  R(T®i  (p))  ^  0  means  that  only  the  contour  points 
of  R  are  involved. 

A  large  number  of  interpolations  are  involved  in  the 
alignment  process.  The  accuracy  of  estimation  of  the  ro¬ 
tation  and  translation  parameters  is  directly  related  to  the 


accuracy  of  the  underlying  interpolation  model.  Simple 
approaches  such  as  the  nearest  neighbor  interpolation  are 

plement,  though  they  produce  images  with  noticeable  arti¬ 
facts.  Besides,  as  the  translation  and  rotation  parameters 
should  compensate  for  accuracy  by  having  subvoxel  values, 
this  type  of  interpolation  would  not  be  appropriate.  More 
satisfactory  results  can  be  obtained  by  small-kernel  cubic 
convolution  techniques,  bilinear,  or  convolution-based  in¬ 
terpolation.  According  to  sampling  theory,  optimal  results 
are  obtained  using  sinus  cardinal  interpolation,  but  at  the 
expense  of  a  high  computational  cost.  As  a  compromise,  a 
bilinear  interpolation  technique  has  been  used  in  the  opti¬ 
mization  steps.  At  the  end  of  the  algorithm,  the  alignment 
parameters  are  refined  using  a  sinus  cardial  interpolation 
that  preserves  the  quality  of  the  image  to  be  aligned.  This 
technique  has  proven  to  be  fast  and  efficient. 

III.  Experimental  Results 


To  evaluate  our  method,  we  applied  the  algorithm  to  the 
reconstruction  of  an  artificially  misaligned  3D  human  skull 
volume  (figure  1).  The  slices  of  the  original  256  x  256  x 
140  CT  volume  were  transformed  using  translations  vary¬ 
ing  from  -10  to  +10  pixels  and  rotations  varying  from  -20 
to  +20  degrees.  The  transformations  for  each  slice  were 
random  following  a  uniform  distribution  in  order  not  to 
privilege  any  slice  (figures  1(a)  and  1(b)).  Table  I  presents 
statistics  on  the  alignment  errors.  The  algorithm  revealed 
robust  in  aligning  this  type  of  image  providing  small  reg¬ 
istration  errors.  Figures  1(c)  and  1(d)  present  the  recon¬ 
structed  volume. 


Fig.  1.  Reconstruction  of  a  3D  human  skull  volume  of  140  slices,  (a) 
Multiplanar  view  of  the  volume  before  registration,  (b)  Three- 
dimensional  view  of  the  volume  before  registration,  (c)  Multipla¬ 
nar  view  of  the  volume  after  registration,  (d)  Three-dimensional 
view  of  the  volume  after  registration. 


Alignment  error  statistics 


A  tx 

A  ty 

A0 

median 

0.23 

0.21 

0.26 

maximum 

1.95 

1.94 

1.64 

mean+s.dev 

0.33+0.32 

0.34+0.33 

0.25+0.25 

TABLE  II 

A  set  of  140  slices  of  a  3D  CT  human  skull  volume  were  artificially 
transformed  using  different  rigid  transformation  parameters.  Each 


Alignment  error  statistics 


A  tx 

A  ty 

A0 

median 

2.10 

0.33 

0.07 

maximum 

1.45 

2.02 

2.42 

mean+s.dev 

0.37+0.28 

0.38+0.30 

0.19+0.35 

slice  was  translated  by  0.15  pixels  in  both  directions  and  rotated  by 
0.3  degrees  with  respect  to  its  preceding  slice.  Different  statistics  on 
the  errors  for  the  rigid  transformation  parameters  are  presented. 
Translation  errors  are  expressed  in  pixels  and  rotation  error  in 
degrees. 


TABLE  I 

A  set  of  140  slices  of  a  3D  CT  human  skull  volume  were  artificially 
transformed  using  different  rigid  transformation  parameters.  Each 
slice  was  randomly  transformed  using  translations  varying  from  -10 
to  +10  pixels  and  rotations  varying  from  -20  to  +20  degrees. 
Different  statistics  on  the  errors  for  the  rigid  transformation 
parameters  are  presented.  Translation  errors  are  expressed  in  pixels 
and  rotation  error  in  degrees. 


Moreover,  we  have  uniformly  transformed  140  slices  of 
the  same  3D  volume  by  applying  to  each  slice  Ii  a  transla¬ 
tion  of  tlx  —  tl~x  +0.15  pixels  and  tly  —  t1^1  +0.15  pixels 
and  a  rotation  of  6l  =  +  0.3  degrees.  As  the  volume 

has  140  slices,  the  last  slice  is  translated  by  21  pixels  in 
both  directions  and  rotated  by  42  degrees  with  respect  to 
its  initial  position.  Table  II  presents  the  registration  er¬ 
rors  of  the  method.  It  is  illustrated  that  our  approach  has 
subpixel  mean  and  median  errors.  Also  maximum  errors 
are  slightly  superior  to  1  pixel  and  1  degree  respectively 
showing  the  robustness  of  the  technique. 


Furthermore,  the  algorithm  was  applied  to  the  recon¬ 
struction  of  volumes  (tooth  germs,  biological  tissues)  with 
unknown  ground  truth.  The  method’s  performance  was 
compared  with  the  manual  alignment  accomplished  by  an 
expert  physician.  Figure  2  shows  the  reconstruction  of  a 
tooth  germ  by  an  expert  dentist  (fig.  2(a)  and  2(b))  and 
by  our  method  (fig.  2(c)  and  2(d)).  It  is  illustrated  that 
human  intervention  fails  to  correctly  align  the  slices,  whilst 
our  method  is  efficient  and  can  achieve  alignment  with  high 
accuracy. 

Also,  Figure  3  depicts  a  3D  tissue  containing  a  large 
number  of  vessels.  Figures  3(a)  and  3(b)  show  the  vol¬ 
ume  aligned  by  an  expert  biologist  and  Figures  3(c)  and 
3(d)  the  tissue  after  alignment  by  our  method.  This  vol¬ 
ume  presents  cuts  and  discontinuities  and  the  tissues  had 
been  stretched  during  the  cut  procedure.  Despite  these 
drawbacks,  according  to  the  expert  biologist,  the  algorithm 
aligned  correctly  the  slices. 

Finally,  let  us  notice  that  the  algorithm  has  a  computa- 


process.  By  these  means,  the  major  problems  set  by  the 
registration  of  serially  acquired  slices  are  addressed.  With 
S  Oct25-2§h^(j)(Jjr) 0f  the  registration  prob¬ 
lem  (rather  than  a  standard  step  by  step,  sequential  formu¬ 
lation)  ,  no  global  offset  nor  error  propagations  are  observed 
in  the  final  alignment.  The  approach  seems  promising  and 
its  association  to  more  sophisticated  but  time  consuming 
_  k  pixel  similarity  metrics  (mutual  information  [18],  robust 
estimation-based  measures  [19])  may  improve  its  accuracy 
and  is  a  perspective  of  this  work. 
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